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flj ■ Abstract. A fundamental hypothesis for the interpretation of the measured large-scale line-of-sight peculiar velocities 

P^ I of galaxies is that the large-scale cosmic flows are irrotational. In order to assess the validity of this assumption, we 
estimate, within the frame of the gravitational instability scenario, the amount of vorticity generated after the first 
shell crossings in large-scale caustics. In the Zel'dovich approximation the first emerging singularities form sheet like 

T-H ■ structures. Here we compute the expectation profile of an initial overdensity under the constraint that it goes through 

^ ' its first shell crossing at the present time. We find that this profile corresponds to rather oblate structures in Lagrangian 

^1^ . space. Assuming the Zel'dovich approximation is still adequate not only at the first stages of the evolution but also 

^, ■ slightly after the first shell crossing, we calculate the size and shape of those caustics and their vorticity content as a 

fvq function of time and for different cosmologies. 

f^ The average vorticity created in these caustics is small: of the order of one (in units of the Hubble constant). 

On ■ To illustrate this point we compute the contribution of such caustics to the probability distribution function of the 

^^ I filtered vorticity at large scales. We find that this contribution that this yields a negligible contribution at the 10 to 

,^ . 15 h~^Mpc scales. It becomes significant only at the scales of 3 to 4 h^^Mpc, that is, slightly above the galaxy cluster 

Q-i" scales. 
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1. Introduction 

The analysis of large-scale cosmic flows has become a very active field in cosmology (see Dekel 1994 for a recent review 
on the subject). The main reason is that it can in principle give access to direct dynamical measurements of various 
quantities of cosmological interest. There are now a very large number of methods and results for the comparison of 
the measured large-scale flows with the measured density fluctuations as observed in the galaxy catalogues. Most of 
these methods are sensitive to a combination of the density of the universe in units of the critical density, fl, and 
the linear bias, b, associated to the mass tracers adopted to estimate the density fluctuations. The estimated values 
of /3 = ri°-^/6 are about 0.3 to 1 depending on the method or on the tracers that are used. There are other lines of 
activities that aim to estimate J7 from only the intrinsic properties of the velocity field, (i.e., without comparison with 
the observed galaxy density fluctuations). All these methods exploit non-Gaussian features expected to appear in the 
velocity fleld, either the maximum expansion rate of the voids (Dekel & Rees 1994), non-Gaussian general features 
as expected from the Zel'dovich approximation (Nusser & Dekel 1993), or the skewness of the velocity divergence 
distribution (Bernardeau et al. 1995). Yet they all also assume that the velocity field is potential. This is indeed a 
necessary requirement for building the whole 3D velocity out of the line-of-sight informations in reconstruction schemes 
such as Potent (Bertschinger & al. 1990, Dekel et al. 1994). This is also a required assumption for carrying calculations 
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in the framework of perturbation theory. It is therefore interesting to check the rotational content of the cosmic flows 
at scales at which they are considered in galaxy catalogues, that is at about 10 to 15ft-^^Mpc. This investigation ought 
to be carried in the frame of the gravitational instability scenario with Gaussian initial conditions. It is known that 
in the single stream regime, primordial vorticity is diluted by the expansion and that the higher order terms in a 
perturbation expansion cannot create "new" vorticity. Hence it is natural to assume that the vorticity on larger scales 
originate from the (rare) regions where multi-streaming occurs. During the formation of large scale structures this 
happens first when the largest caustics cross the first singularity, creating a three-flow region where vorticity can be 
generated. As we argue in Sect. 2, analytical calculations of constrained random Gaussian flelds suggest that the largest 
caustics that are created are sheet-like structures, in rough agreement with what is found in numerical simulations or 
in galaxy catalogues. It is therefore reasonable to use Zel'dovich's approximation to describe the subsequent evolution 
of those objects. 

In order to estimate the large scales vorticity distribution we therefore proceed in flve steps: flrst we evaluate the 
mean constrained random field corresponding to a local asymmetry of the deformation tensor on a given scale, Rl] 
secondly we solve for the multi-flow regime within the generated caustic, using Zel'dovich's approximation throughout, 
even slightly beyond this first singularity. We then evaluate the vorticity field in that caustic. The next step involves 
modelling the variation of the characteristics of typical caustics as a function of time for different power spectra. 
Finally, we estimate the amount of vorticity expected at large scales arising from large scale flow caustics. 

For the sake of simplicity and because is pedagologically more appealing, we present calculations carried out in 
two dimensions as well as in three dimensions. The former case is in particular easier to handle numerically. 

The second Section of this paper evaluates the characteristics of the typical caustics expected at large-scale in a 
2D or 3D density field. The third Section is devoted to the explicit calculation of the vorticity for the most typical 
caustics. The fourth Section provides an estimate for the shape of the tail of the probability distribution function of 
the modulus of the vorticity in a sphere of a given radius. It is followed by a discussion on the validity and implications 
of these results. 

2. Asymmetric constrained random fields 

Since it is not our ambition to solve the problem of deriving the vorticity statistics in its whole generality the vorticity 
will be estimated only within specific but typical caustics in the framework of the Zel'dovich approximation. 

The first step involves building an initial density field in which a caustic will eventually appear. The initial fluctu- 
ations are assumed to be Gaussian with a given power spectrum P(fc), characterizing the amplitude and shape of the 
initial fluctuations. No a priori assumptions about the values of H, and A are made. It will be shown that the statistics 
has very straightforward dependences upon these parameters. The expectation values of the random variables, (5(k), 
corresponding to the Fourier transforms of the local density field, 

(5(x) = /"d3k(5(k) exp[ik-x], (1) 

are calculated once a local constraint has been imposed. This constraint will be chosen so that the caustic-to-be will 
have just gone through first shell crossing at the present time. It is expressed in terms of the local deformation matrix 
of the smoothed density field. The components of the local deformation tensor at the position Xq are given by 

$«.,(xo) = /"d3k(5(k) VKz5(fci?L)exp[ik.xo]^, (2) 

where Wd is the adopted window function. In what follows, we will use the top- hat window function for which 

^2(fc) = 2^ in 2D, W^{k)=i./^^-^0^ in iD, (3) 

where J,y is the Bessel function of index v. The scale Rl is the scale of the caustic in Lagrangian space. Here ctq stands 
for the rms density fluctuation at this scale: 

al^ Jd'kP{k)Wi{kRL). (4) 

For the sake of simplicity a typical caustic is chosen to be characterized by the average local perturbation over a 
sphere of radius Rl for which the deformation tensor at its centre is fixed. We are aware that this is a somewhat 
drastic approximation but consider that, at large scales, the behaviour of caustics having the mean initial profile will 
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be representative of the average behaviour. This is certainly not true at small scales where the conrplex interactions of 
structures at different scales and positions are expected to affect the global behaviour of any given caustic. For some 
rare enough objects however we expect the fluctuations around the mean profile to be small enough to affect only 
weakly the global properties of the caustics. This has been shown to be true in the early stages of the dynamics for 
spherically symmetric perturbations (Bernardeau 1994a). In the following we will, however, encounter properties (see 
§3.3) that we think are not robust against small scale fluctuations. Such properties will be ignored in the subsequent 
applications of our results. 

Within the frame of this calculation, the values of (5(k) hence correspond to the expectation values of (5(k) for the 
power spectrum P{k) when the constraints on the deformation tensor are satisfied. These solutions can be written as 
a linear combination of the values of the deformation tensor: 
D (C^^) ^ 



(C-^)o,o 



i=l 



where the coefficients C is the matrix of the cross-correlations between the random Gaussian variables <I>y and i5(k) 
as shown in Appendix A. In Eq. (5) the summation is made only on the diagonal elements of the deformation tensor 
since it is always possible to choose the axis in such a way that the other elements are zero. In this instance, the 
diagonal elements are identified with the eigenvalues Ai, of the matrix. 

2.1. The 2D field 

In 2D geometry, the two coefficients ai and a2 defined by Eq. (5) are given by 

ai = (3/i - /2)/a2 , a2 = (3/2 - /i)/(t2 , where h ^ {Sk<^^^) ^ P{k)W2{kRL)'^. (6) 

The brackets, /.V denote ensemble averages over the initial (unconstrained) random density field. As a result, Eq. (5) 
reads 

<5(k) = ^W^2(fcflL) ^2 (^^ + A2) + 4 (Al - A2) cos(2e)] ; (7) 

Al and A2 are the eigenvalues of the deformation tensor and where 6 is the angle between k and the eigenvector 
associated with the first eigenvalue (see Appendix A for details). Consider the parameter a defined by 

Al + A2 

The coefficient a represents the amount of asymmetry in the fluctuation (thus a = corresponds to a spherically 
symmetric perturbation). This parameter is similar to the eccentricity, e, that was used by Bardeen et al. (1986) and 
more specifically by Bond & Efstathiou (1987) for 2D fields. In these studies however investigations were made for the 
shape of the peaks around the maximum (i.e. eigenvalues of the second order derivatives of the local density), so a 
and e cannot be straightforwardly identified. 

The formation time of the first singularity is determined by the maximum value of the eigenvalues, Amax- It is 
therefore relevant to calculate the distribution function of Amax, and the distribution function of a once Amax is known. 
From the statistical properties of the matrix elements ^ij we derive the distribution function of the eigenvalues Amin 
and Amax (see Appendix B), which reads 

1 /3 

with 

Jl = Amin + Amax , ^2 = Amin Amax- (10) 

The distribution function of Amax follows by numerical integration over Amin- Fig. (1) shows the distribution function 
of Amax in units of the variance. The dashed line corresponds to the approximation, valid at Amax/co 3> 1: 

n A 

(11) 



23/2 

-'^l Amin, '^niaxl — 77o T lAniax -^rain) CXp 
TT-^/^ (To 



'0 



;2io^i-4J. 



(9) 



Aniax 



Pmax(Amaxj CiAmax ~ 1-5 CXp 



4 

3 V o-Q 



CTO 



The distribution function of a for different values of Amax/co is presented in Fig. (2). It turns out that the most 
significant value corresponds to a w 1. In the following this value is chosen as the typical value for the asymmetry in 
two dimensions. 
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Fig. 2. The distribution functions of a for fixed values of 
Fig. 1. The distribution function of Amax/o"o (solid line) Amax/o"o = 1, 2, 3, 4 (respectively the solid, long dashed, 
in 2D dynamics. The dashed line is given by (Eq. (11)): short dashed and long dotted dashed lines). 



2.2. The 3D field 

In three dimensions the geometry is slightly more comphcated and yields for the constrained density field (see Ap- 
pendix B for details) 



<5(k) 



iP{k)Wz{kRL) 



%al 



(\i [1 + 5 cos(20fe) - 5cos(26'fc) - 5 cos(2(/)fe) cos(20fe)] 



Aa [1 + 5 cos(20fe) - 5 cos(26'fe) - 5 cos(20fe) 008(26*^:)] + 2A3 [3 + 5 cos(2 Ok)] ) , 



(12) 



where Ok and 0^ are polar angles of the vector k with respect to the basis of the eigenvectors associated to the three 
eigenvalues, Ai, A2, A3. The asymmetry of the distribution is again characterized by the values of 



2A3 — Al — A2 

a = 5 — , and = 5 ■ 

Al + A2 + 6A3 



Al — A2 



Al + A2 + 6A3 



(13) 



When b only is zero Eq. (13) corresponds to a perturbation with axial symmetry, and when both a and b are zero it 
is a spherically symmetric perturbation. In terms of a and b Eq. (12) then becomes 



<5(k) = 3^(fc)W'3(fcJ?L) ^^^ ^ ^^ ^ g^^^ A ^ ^ cos(20fc) + b cos{2cj)k) [l + cos{29k)] 



(14) 



Let us now evaluate the distribution of a and b from the distribution function of the eigenvalues (Ai,A2, A3) in 3D 
(assuming Ai > A2 > A3) in order to identify the shape of the most significant caustics. This distribution is given by 
(Doroshkevich 1970) 



5'''/2 27 
P(Ai, A2, A3) = ^^^ (Al - A2) (Al - A3) (A2 - A3) cxp 



Sttct" 



3Jt 



^ '^J. 



(15) 



with 



Ji = Al + A2 + A3 , and J2 = A1A2 + A2A3 + A3A1. 



(16) 



From this expression we compute numerically the distribution function of the maximum eigenvalue (Fig. (3)). An 
analytical fit of this distribution function is provided by its behaviour at large Amax 



Pma^y^maKJ *^^n 



6 



A. 



c^o 



exp 



5 /X 

'2 



CTO 



dAn 



CTO 



(17) 
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This fit is accurate for the rare event tail (as shown in Fig. (3)), which wiU be relevant for the derivation of Sect. 4.4. 
For a given value of Amax we compute the distribution of the other eigenvalues, and thus the join distribution function 
of a and b. 
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Fig. 4. The contour plot for the distribution of a and 
b for a fixed value of Amax/o"o ~ 2 (dashed lines) and 

Fig. 3. The distribution function of Amax/cro (solid line) Amax/o"o — 3 (solid lines). The lines are evenly distributed 

in 3D dynamics. The dashed line is the analytical fit (17). in a logarithmic scale. 

The resulting contour plot corresponding to Amax/co — 2 and Amax/co = 3 is illustrated on Fig. (4). As for the 
distribution of a in the previous subsection in 2D it depends only weakly upon the adopted value of Amax (although 
the position of the maximum varies a little), and it tends to be all the more peaked on its maximum as Amax is large. 
This implies that a typical caustic will be given by a f* 1 with a small b. For further simplifications we will assume that 
6 = 0. Such a caustic then corresponds to a pancake-like structure with axial symmetry. Note that this result seems 
to differ from the results of Bardeen et al. (1986) who found that the shape of the rare peaks should be somewhat 
spherically symmetric or filamentary (this picture was recently sustained by Pogosyan & al., 1996, from the result of 
iV-body simulations). This apparent discrepancy is due to the constraint under which the expectation values of a and 
b are calculated. In Bardeen et al.'s work the constraint is given by the value of the local density, i.e. the sum of the 
three eigenvalues, whereas in this paper we put a constraint on the largest eigenvalue. This is a natural assumption 
for this investigation since the multi-streaming occurs as soon as a singularity has been reached in one direction. Of 
course, this analysis assumes that the Zel'dovich approximation holds in order to predict the time at which this first 
singularity is reached. For oblate initial structures such as the ones obtained for the most likely values of a (see Figs. 
5 and 6), we expect that this approximation is sufficiently accurate. 

3. The geometry & vorticity of large-scale caustics 

In this section we investigate the properties of the caustics that are induced by the initial density fiuctuation profiles we 
found in the previous section. All the calculations are performed within the framework of the Zel'dovich approximation, 
even sightly after the first shell crossing. 



3.1. The linear displacement field 

In the framework of the Zel'dovich approximation the displacement field can be written 

^ = q + D{t)/D{to)^icD; 



(18) 



where D{t) accounts for the time dependency of the linear growing mode (it is proportional to the expansion factor 
in case of an Einstein-de Sitter geometry only). An important simplification is that, at the order of the Zel'dovich 
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approximation, this displacement field is separable in time and space, and its space dependence, ^(q), is potential, 
i.e., there is a velocity potential U{q) so that 

*(q) = Vq • ;7(q) . (19) 

This velocity potential is given by 

C/(q) = y d3q<5(k) -1 exppk • (q - qo)] . (20) 

By construction the point qo in Lagrangian space corresponds to the point xq in real space (central position of the 
caustic). Both of them will be taken to be zero. For the calculation of the explicit expressions of (5(k) and U{q) we 
will assume that the power spectrum follows a power law behaviour, 

P{k) ex fc", (21) 

characterized by the power index n. From Eq. (21) the expression of the linear variance as a function of scale follows 

a{RL) oc i?-("+^)/^ (22) 

This approximation is valid within a limited scale range as will be discussed in Sect. 5. At the scales of interest the 
index n is expected to be the range of n « — 1, —2 from the constraints obtained with the large-scale galaxy catalogues, 
like the APM survey (Peacock 1991) the IRAS galaxy survey (Fisher et al. 1993) or from X-ray cluster number counts 
(Henry & Arnaud 1991, Eke et al. 1996, Oukbir & Blanchard 1997). In two dimensions there are of course no such 
observationally motivated values, but we will consider n of the order of —1 as an illustrative case. 

3.1.1. The 2D potential 

From the Eqs. (7), (20) it is possible to calculate the expression of the potential 

U{q) = G(0,n-2,q)+acos(29„)\G(0,n-2,q)-2G{l,n-2,q)] , with G(iy,n,q)= /d^kfc" ^7^ W^2D(fc)-(23) 

J [kqr 

The latter expression is given by 

Giiy,n,q)= 2Fi{l + n/2,n/2,l + u,q^), for q<l, and (24) 

G(^, n, q) = ^^it2 w^^^ ~ %^^ 2^1 (1 + n/2, l-v + n/2, 2, q-^) , for q > 1. (25) 

The expressions for the gradients of the potential involve similar hyper-geometric functions. 

3.1.2. The 3D potential 

The expression of the potential following from Eqs. (12), (20) becomes quite complicated, but involves here only 
"simple" functions. It reads 

C/(q) = [V{q) - V{-q)]/q\ (26) 

with 

V{q) = \l + q\'^-'^s\gn{l + q)[A{q)-B{q)[hcos{2(t)) [1 - cos(2 g)] + a cos(2 61)] ) , (27) 

A{q) = -IQq^ + 7nq^ -n^ q^ + ^q^ -nq^ + a {-l + 2q-nq + 2q^ -nq^ - q^) , (28) 

B{q) = S-Gq + Suq + Aq"^ -Anq"^ + n'^q^ -2q^ + nq^ (29) 

Note that the potentials in Eqs. (23) and (26) have discontinuous derivative at g = 1, which is an artifact of using 
a top-hat window function. Note also that the potentials given here have arbitrary normalizations. This is of no 
consequence for the derived results since the global normalization of the initial density profile is absorbed in the 
discussion for the value of Amax (Sect. 4.4). 
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Fig. 5. The shape of the caustic for the 2D dynamics, Fig. 6. The shape of the caustic for the 3D dynamics, 
n = —1, and Amax ~ 1-3. The dashed line is the shape n = —1.5 and Amax ~ 1-5. The external shell is the La- 
in Lagrangian space and the solid line the shape in real grangian position of the caustic, the internal one its posi- 
space. tion in real space. 



3. 2. The shape of the caustics 

A multi-flow region forms as soon as Eq. (18) has more than one solution. The corresponding region forms the so-called 
caustic. These regions are illustrated in Figs. (5) and (6) in respectively 2 and 3 dimensions for typical values of the 
parameters. The solid lines show in 2D the shape of the caustic in real space, and the dashed lines their shape in the 
original Lagrangian space. 

For the chosen values of a and b and for the relevant Amax the caustics form elongated structures. These figures 
are plotted in units of the smoothing scale Rl ■ They suggest that the largest dimension of the caustics are roughly of 
the order of magnitude of the initial Lagrangian scale. Note that the boundaries of the caustics correspond to surfaces 
(or lines in 2D) where the Jacobian of the transformation between Lagrangian space and real space vanishes, i.e. 



^(q) 



I — I 
' 9q' 



= 0. 



(30) 



The size and shape of these caustics are characterized, in 2D and 3D (although only approximately), by two lengths, 
the half-depth of the caustic, d, (that is the distance that has been covered by the shock front after the first singularity) 
and its half-extension e. For instance in Fig. (5) the value of d is about 0.1 and the value of e is about 0.9 in units of 
the Lagrangian size of the fiuctuation R^ ■ In the case of the 3D dynamics e corresponds to the radius of the caustic 
since we restrict ourselves to cylindrical symmetry. 

The density in each fiow "s" is given by the inverse of the Jacobian of the transformation so that 



p(q,) = l/J(q,). 

The total density within the caustic is then given by the summation over each flow of each of their densities, 



(31) 



(32) 



flow s 



3.3. The velocity field, and the generated vorticity 
The velocity in each flow is given by 

u{q) = D{t)/D{to)^{q). 



(33) 
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Fig. 7. The map of the vorticity in a typical 2D caustic {n — —1). Left panel: the local vorticity is antisymmetric with respect 
to the centre of the caustic. It points along the Z-axis, and is positive in the second and fourth quadrant, and negative in the 
first and third. Right panels: behavior of the local vorticity along two different lines (thick dot-dashed line on the left panel). 
The top panel shows that the vorticity is singular near the edge of the caustic. It behaves as described by Eq. (43) and there is 
a non zero lineic vorticity located on the edges (represented here by a vertical line) due to the discontinuity of the local velocity 
field. The bottom panel shows that the local vorticity goes continuously to zero towards the axes. 



For a given Robertson Walker cosmology, D{t) obeys 
D{t) = f{n) Ho D{t) « f^of' Ha D{t) . 



(34) 



where Hq is the Hubble constant at the present time and /(fi) is the logarithmic derivative of the growing factor with 
respect to the expansion factor. Eq. (34) is the only place where the Vt dependence (and A dependence though it is 
negligible) will come into play. 

In general the velocity field, u(x), is defined as the density averaged velocities of each flow. Thus we have 



u(x) 



EflowsP(q^)"(q«) 
Eflow.p(q^) 



(35) 



where the summation is carried on all the flows that have entered the neighborhood of x. The vorticity is then given 
by the anti-symmetric derivatives of the total velocity with respect to x: 



C..(X) . ^,'=..^=5^,'=.. IW ^(i^-),,U.(q.)l W p(q.) 



^ p(q^)u,(q,) 



.flow s 



E 

.flow s 






I 



.flow s 



where Di^^ is the matrix of the transformation between the Lagrangian space and the Enlerian space, 
D . -^ 

aqj 



(36) 



(37) 



and (^-'^'"^ the totally antisymmetric tensor. The numerical expression of the local vorticity follows from the roots of 
Eq. (18) and the potentials Eqs. (23), (26). 
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3.3.1. The local vorticity 



As illustrated in Fig. (7) (the 2D case) and (8) (the 3D case), the vorticity is null outside the caustic. First note that 
the vorticity sign changes from one quadrant to another, so that the global vorticity is zero (as it should be), and note 
that within each quadrant the vorticity is rather smooth. Note also that the vorticity is mainly located near the edges 
of the caustic. In fact the vorticity at the edge is unbounded and the behaviour of the vorticity close to the edges is 
easily estimated. Calling qo and Xq the position of a point on the edge in respectively the Lagrangian space and the 
Eulerian space, we can expand x and q close to Xq and qo. Since the linear term in the expansion is singular in q = qo 
(by definition of the caustic), there is one direction, orthogonal to the edge and typeset with the subscript _l, for which 

(x-Xo)_L w-?7(qi-qo)i, (38) 

where ry is given by the second order expansion of the displacement field along this direction. The minus sign accounts 
here for the fact that xo_l has been assumed to be larger than xj^. This equation is valid for two different flows (say 1 
and 2) corresponding to the two roots of q^ in Eq. (38). The Jacobian for the first two flows is then 



J(x) w -277 (q« - qo)± ~ 2 v/?7(xo-x)i. (39) 

Note that on the edge of the caustic, J(x)|9 J(x)/9x| has a finite value, r/. There is also a third flow in the vicinity of 
xo which is regular; let us call qa the Lagrangian position of Xq in this flow. The velocity is then given by 



u(x) « (^(xo - x)-^'7v^u(qo) + p(q3) u(q3) j / (^(xq - ^)l'^VVv + ^(qa) j ■ (40) 

As a result we have 

u(x) w u(qo) + p(q3) V^ (xo - x)^^ (u(q3) - u(qo)), (41) 

when x is within the caustic and 

u(x) = u(q3), (42) 

when X has crossed the caustic boundary. The local velocity is thus discontinuous at the caustic boundary and the 
induced vorticity is consequently singular at xo with 

^(x) « -p(q3) V^(xo - x)-^/' (u(q3) - u(qo))||/2. (43) 

The direction || is a direction parallel to the caustic. There is only one such direction in 2D, two in 3D. There is 
however not only a surface (or volume) contribution within the caustic. Because of the discontinuity of the velocity 
fleld at the edges of the caustic, a vorticity field on the boundary of the caustic is created (see Fig. 7 for the 2D case), 
whose linear or surface density for respectively the 2D and 3D cases are given by 

t^lin., surf = (u(q3) -u(qo))||. (44) 

It turns out that the two contributions tend to cancel each other. Indeed, as we have noticed previously, the velocity 
increases close to the edge of the caustic, and then has a discontinuity at the edge. This creates a sharp peak in the 
vicinity of the edge of the vorticity. The vorticity, which is obtained by differentiation of the local velocity is then 
expected to be opposite on both side of this peak. Realistically, the small scale perturbations are going to wash out 
these features, and to smooth the velocity peaks. As a result the quantities describing the behaviour of the vorticity 
near the edge of the caustic are not robust and should not be taken at face value. On the other hand, we expect the 
integrated vorticity to be a more robust quantity, since it is roughly independent of small scale fluctuations. 

3.3.2. The integrated vorticity 

In two dimensions, the integrated vorticity in each quadrant can be easily obtained numerically by simple one dimen- 
sional integrals which, from Stoke's theorem, can be expressed as 

t^quad. = / d^xw(x) = / u-dl, (45) 

Jquadran J edges 
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where dl describes the edge of the quadrant. One should bear in mind that, in Eq. (45) the velocities on the edge of 
the caustic are taken as the velocities of the third flow, u(q3), so that the singular part of the vorticity is taken into 
account. 

In three dimensions and for (almost) spherically symmetric caustics the local vorticity is independent of the 
azimuthal angle, 9. It is then convenient to calculate the integrated vorticity per azimuthal angle in each quadrant. 



w, 



quad. 



de 



dzrdrw(x) I dO 



quadran 



•udl 



edges 



d xw. 



quadran 



AO, 



(46) 



where r is the distance of the running point to the symmetry axis, and Uz is the velocity component along this axis. 
Compared to the 2D case there is a further difficulty due to the surface integral of one component of the velocity. 
Note nonetheless that this contribution is not singular at the edge of the caustic as shown by Eq. (41), and can thus 
be safely computed numerically. We found that this second integral contributes typically to about 15% of the total for 
the relevant caustics. 



0.10 

: Height 

0.05^ 

0.00^ 

-0.05 ^ 

-0.10^ 
-1.0 




-0.5 Radius 



1.0 



Fig. 8. Section of the vorticity field for the caustic of Fig. (6). The local vorticity is antisymmetric with respect to the centre 
of the caustic. In this X — Z section, it points along the Y-axis, and is negative in the second and fourth quadrant, positive in 
the first and third. 



3.3.3. Scaling laws 

We now bring forward fits to describe the dependence of the integrated vorticity with the spectral index n and Amax. 
which will allow us to characterize the most significant caustics that contribute to the large-scales vorticity. We make 
explicit the dependence of those quantities with respect to the size of the perturbation R^ and the cosmological param- 
eter J7. Expressed in units of the expansion factor, the displacement, in the Zel'dovich approximation, is independent 
of Vt. Therefore a and h are independent of fi, and are simply proportional to R^. The total vorticity in each quadrant 
is on the other hand proportional to Hq and /(SI) (defined in Eq. (34)), given that it is proportional to the local 
velocity, and is clearly proportional to the volume of the perturbation. We thus have the following scalings. 



d{RL) = Rl do (A„ 



l)"^ e(i?L) = Rl eo (A,nax - 1)"% u;,^,,,d{RL,^) 



fin) i?f C^o (Amax - 1)' 



Hn 



(47) 



where the parameters a, a^i, ae, ujq, do and eo are given in Table (1) and (2) for respectively the 2D and the 3D 
geometry. The accuracy of these fits is illustrated on Figs. (9)- (10). These functions yield estimates of the geometry 
and vorticity generated by these large-scale caustics. From these tables one can see that the average vorticity (in units 
of Hq) is roughly one within the caustic. The amount of vorticity which is generated in the caustics is thus found to 
be somewhat limited. It is also interesting to note that LUquad. presents no singular behaviour when the caustic appears 



at An 
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I.e. 



a>l). 



4. The vorticity distribution at large scales 

As argued previously, the calculation of the global shape of the vorticity distribution is beyond the scope of this paper. 
Indeed the low lo behaviour of the vorticity distribution is dominated by the small caustics that are not rare, and 
therefore not well described by the dynamical evolution of an isolated object. The aim of this section is to estimate 
the shape and position of the cutoff in the probabihty distribution function of the local smoothed vorticity. We will 
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Fig. 9. ci;quad. for 2D caustics as a function of Amax and Fig. 10. cjquad. for 3D caustics as a function of Amax and 

its corresponding fit for a n = —1.5 (circles, solid line), its corresponding fit for a n = —2 (circles, solid line), 

n = —1 (squares, long dash line), and n = —0.5 (triangles, n = —1.5 (squares, long dash line), and n = —1 (triangles, 

short dashed line) power spectrum. short dashed line) power spectrum. 

Table 1. Fitting parameters in Eq. (47) for the 2D caustics. The quality of those fits for loq and a are illustrated in Fig. (9). 
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Table 2. Fitting parameters in Eq. (47) for the 3D caustics. The quality of those fits are illustrated in Fig. (10). 
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therefore estimate Pr^ (> w,,), the probability that in a circular or spherical cell of radius Rg the mean vorticity exceeds 
ujs- This estimation requires 

(i) identifying the caustics that contribute mostly for each case; 
(ii) estimating the contribution of each of those caustics. 

In each case various approximations are used. In the main text we simply spell the major highlights of the derivation. 
A more detailed and explicit calculation of the vorticity distribution is presented in Appendix C. 



4-.1. Identification of the caustics 

We assume in what follows that LOg is large enough for the contribution to Pr^{> uJs) to be dominated by large and 
rare caustics. This assumption is the corner stone of the calculation: only a small fraction of the caustics with specific 
characteristics at some critical time will contribute. 

The identification of the caustics contributing most results of a trade off between the amount of vorticity a given 
caustic can generate and its relative rarity: the higher Amaxi the greater the internal vorticity is, according to Eq. (47) 
and given that a is positive, but the rarer those caustics are (Equations (11) and (17)). Obviously Amax should be 
larger than unity for any vorticity at all to be generated. The calculation is slightly complicated by the fact that the 
Eulerian size of the caustics also depends of the value of Amax- Let us assume here that the Eulerian size of the caustics 
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is substantially smaller than the smoothing length, so that the entire integrated vorticity in a quadrant can contribute 
(in Appendix C, this assumption is shown to be self-consistent). This implies a scaling relation between the smoothing 
cell, Us and Amax, 

UJs Rf (X i?f (A„,ax - 1)". (48) 

For a given smoothing length and a given ujs, Eq. (48) yields a relation between the value of Amax and the size of the 
caustic. The caustics which contribute most to the vorticity are then obtained by minimizing the ratio X^^^/a'^{RL) 
which appears in the exponential cutoff of the distribution function of Amax (Equations (11) and (17)). Given that 
<7^(Rl) behaves like -R^ this minimization yields for the extremum value of Amax, 

A(") = — (49) 

■""" 2D -a{n + D)' ^ ' 

Note that for the values of a we have found, Amax is always finite and positive. This means that the filtered vorticity 
is indeed expected to be dominated by caustics which have evolved for a finite time. This provides an a posteriori 
justification of the assumptions leading to this calculation. 

The value of Amax found in Eq. (49) is a robust result of our calculations, although it cannot be excluded that this 
value could be affected by the failure of the Zel'dovich approximation after the first shell crossing. 

4-. 2. Estimation of the caustic contribution to the vorticity PDF 

In order to estimate the contribution of those caustics to Pr^{> Ug) two other fundamental quantities have to be 
estimated: 

(i) the number density of caustics; 
(ii) the volume for which each of them contributes to Pr^ (> lUs). 

These quantities have been estimated for the specific caustics we have previously identified in Sect. 4 1. 

4.2.1. The number density of caustics 

Estimating the number density of caustics is, in general, a complicated problem. In the case of Gaussian fields the 
corresponding investigation was carried by Bardeen et al. (1986) for 3D fields, and by Bond & Efstathiou (1987) for 2D 
fields. The number of caustics is simply determined by the number of points at which the first derivatives of the local 
density vanishes. This defines accordingly the extrema of the local density field. The further requirements we have here 
on the second order derivatives of the potential ensures that such points are in fact maxima of density field. We refer 
here to Bardeen et al. (1986) for more details on how to carry the investigation. A critical step involves transforming 
the (^Dirac function in the value of the first derivatives into a (^Dirac function in the position, thus introducing the 
Jacobian of the second order derivatives of the density field. After some algebra we find: 

n«J{AJ)d A.=p(^|-^^|j ^z,(^^) (2,,2)z,/2 ' (50) 

where the probability p is given either by Eq. (9) or (15) in respectively 2D and 3D, Jac2({Ai}) is the Jacobian of the 
second order derivatives of the density field for given eigenvalues of the deformation matrix and cti is the variance of 
the derivatives of the local density field, 

al{RL) = /d^fc P{k) y WURl)- (51) 

For a given geometry {i.e. given values of a and h) Jac2 is proportional to Am^^,^, and it scales as R^ due to the 
derivatives involved in the expression of the matrix elements. It is therefore appropriate to re-express Eq. (50) as 



n\^\AU^ I) ^^ W d^-^» no({Ai}) /Amax\^ , ,.,.. |Jac2({AJ)| '"-"'^ 

n^J{A.})d A.^p — — T71?-T ^here no({A.}) = 



a(i?i)/; a^(i?i) i?f WRl)) "''''' ASa. (2^)^/2 



CTi 



i?f.(52) 



Note that no, thanks to the prefactor i?£, is a dimensionless quantity in Eq. (52). A further simplification is provided 
by the fact that for large enough values of Amax, the distribution function p({Ai}), at fixed Amax, allows only a small 
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range of possible values for the smaller eigenvalues. We therefore neglect the variations of Jac2({Ai}) with respect 
to those variables: it is viewed here as a function of Amax only and calculated for fixed values of the a-symmetry 
parameters a and b. The ratio cr/ai depends only on the value of the power law index. Recall however (see Bardeen 
et al. 1986) that this ratio is not well-defined for top-hat window functions because of spurious divergences for some 
values of n. To avoid this problem, we used the Gaussian window function to compute this ratio. As a result, for 
fixed values of a and b, uq is a dimensionless quantity that can be explicitly calculated in a straightforward manner. 
Relevant values of rig are given in tables in the Appendix C. 

4-.3. The contributing region 

The region over which each caustic contributes is the surface (or volume in 3D) of space in the vicinity of a given 
caustic where, if one centers a cell in that location, the total vorticity induced by the caustic within the cell is above 

LOs. 

In general the contributing surface or volume can be written, 
Kaus.(i?L, i?«, {AJ, u;,)^ Je K (c, Rl, Rs, {A.}) - w,] d^c, (53) 

where & is the Heaviside step function, c stands for the vector pointing to the center of the sampling sphere, while loc 
is the vorticity found in that sphere intersecting the caustic with characteristics Rs,{Xi}. In its full generality, t4aus. 
is a rather complex function of the scales Rl and Rg , and the eigenvalues Ai through the shape of the caustics and 
of LOs- Yet, since the functional form of the rare event tail in the probability distribution function is basically fixed by 
the exponential in Eq. (11), the only required ingredient for computing Pr^{> lOs) is the scaling behaviour of K;aus. 
at its takeoff - when reaching the critical time, Amax, at which a given caustic is large enough to start contributing. 
The detailed geometry of the caustic and its vorticity field accounts only for a correction in a multiplicative factor. 
Consequently we make approximations describing the distribution of the vorticity on the caustic in order to estimate 
the scaling properties of V^aus. • 

4.3.1. The 2D contributing surface 

In two dimensions we make the radical assumption that the vorticity is entirely concentrated on four discrete points, 
which - consistently with the hypothesis of Sect. 3.3.2, have been taken to bear either the vorticity +a;quad. or — Wquad.j 
depending on which quadrant is being considered. In practice the position of the points is chosen somewhat arbitrarily 
at a third of the depth and extension of the caustic. The corresponding area Vcaus. is therefore identically null before a 
critical time corresponding to the chosen w<, and Rs and then takes a constant value which can be deduced geometrically 
from the area of the loci of the center of the sampling disks. In Fig. (11) we show the shape of this location on a 
particular example. Under this assumption, the function T^aus. takes the form, 

l^caus. = Vo{Rl/Rs) e(Amax - A^L) Rl Rs , (54) 

where Vq can be calculated for the values of interest of Rl and Rg ■ 

4.3.2. The 3D contributing volume 

In three dimensions, the vorticity will be assumed to be distributed uniformly along two rings which are taken to bear 
the linear vorticity 3cijquad./e - with respectively prograde and retrograde orientation. In practice these rings are also 
positioned at a third of the depth and extension of the caustic. The mean vorticity to be expected in a sampling sphere 
of radius Rg is then given by algebraic summation over the segments corresponding to the intersection of that sphere 
with the two rings. Maps of the sampled vorticity as a function of the centers of the sphere are derived to compute 
K;aus. which according to Eq. (53) corresponds to the volume in space defined by these centers yielding a vorticity 
larger than tOs- Fig. (12) gives the shape of this location for a given caustic and sampling radius. The function Vcaus. 
takes the form, 

X^caus. = Vo{Rl/Rs)RlR^s (Amax - A^L)"' , (55) 

where Vq and 7 can be calculated for the values of interest of R^ and Rs at this critical values (see Appendix D, 
where it is in particular demonstrated that when Rl <C Rs, Vq asymptotes to a fixed value and 7 = 1/2). 
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Fig. 11. Sketch showing the adopted simplification for 
describing a 2D caustic. Vorticity is assumed to be local- 
ized on the black dots having either -t-fi^quad. or — ciJquad.- 
The dashed area represents Vcaus. for ujs > jti^quad.l- 



Fig. 12. Sketch showing the adopted simplification for 
describing a 3D caustic. Vorticity is assumed to be lo- 
calized on two rings (that appear as two horizontal 
black segments) having a lineic vorticity of either either 
+3ciJquad./s or — Sojquad./fi- The shaded area represents 
dVci,uB./dujs- 



4-. 4- Estimation of Pj^^{> u>s) 

The tail of the probability distribution for the vorticity is now estimated while integrating over all the caustics that 
might contribute, and assuming that, for a fixed caustic, the probability distribution is given by the number density 
of caustics times the volume associated with each caustic. There is however a further difficulty. The distribution of 
caustics UR^ is well defined for a fixed value of Rl only, but there are actually caustics of all sizes. To circumvent this 
difficulty we simply choose Rl so that the result we obtain is maximal, i.e.. 



PRsi> ^s) — max 

Rl 



d^A, nR,{{X,}) yeaus.(i?L,i?.,{A.},Ws) 



(56) 



Furthermore, it is fair to neglect the dependence of no(Ai) and 14aus on the initial asymmetry because the overall 
factor p{Xi) peaks in a narrow range of relevant values for the smaller eigenvalue(s). It is then possible to integrate 
over those variables introducing the probability distribution of A^ax in the expression of Pr^{> uJs), 
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Pniax l^-^niax J 



no(A„ 
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ctIRl) 
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max? ^^s ^ 



(57) 



We show in Appendix C that the maximum of Eq. (56) is indeed given by caustics of size of the order of Rs at 
most. A detailed account of how to perform the sum in Eq. (56) is also given there for the two geometries. Repeated 
use of the rare event approximation together with the geometrical assumptions on the vorticity distribution sketched 
in Sect. 4.3.1 and Sect. 4.3.2 yields eventually an explicit expression for the tail of the probability distribution for the 
vorticity as a function of ujs and Rg- 



4.4.1. The two dimensional vorticity distribution 

In two dimensions, the vorticity distribution is shown to obey (Eq. (C9)) 



PR^{>iUs)^OMnoVo 



X 



(0) 
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(58) 
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In the rare event regime, the quantity that dominates Eq. (58) arises from the exponential cutoff. For n = — 1 wc find 
for instance that 

1/2 

\og[Pn^{>u;s)]o,3.5^f—. (59) 

The r.h.s. of Eq. (59) is roughly 0.5 when utg « 10^'^, <j{Rs) ~ 0.5 or LOg ~ 0.1, (y{Rs) ~ 1-5, hence defining a threshold 
corresponding to a one sigma damping for Pr^{> Us). Equation (58) is illustrated on Fig. (13). 

4.4.2. The three dimensional vorticity distribution 

Similarly, the probability distribution is shown in the Appendix C (Eq. (C19)) to obey in 3D: 

2 



PRA>^s)^OA8noVo 





jn+3^(n+3)/3 



(60) 



For n = —1.5, Eq. (60) gives for log[Pfl^(> iVg)] 

T2(i?,)- 



log [Pr^ (> uj,)] ~ 20. 4^. (61) 



yielding again at a one sigma level the range of relevant values for ujs and (t(Rs): uJs ~ 510^'"', (t{Rs) ~ 0.5 or 
LUg ~ 0.1, cr(i?s) ~ 3.5. In both cases the caustics start to generate significant vorticity only at rather small scales. 
Equation (60) is also illustrated on Fig. (14). From this figure it is clear that the amount of vorticity that we derived 
is below what has been measured in A^-body simulations (open and filled circles). Numerical measurements of this 
quantity are sparse, so we compared our estimations to measurements carried out by Bernardeau & Van de Weygaert 
(1996) in an adaptive P^M simulation with CDM initial conditions (see Couchman 1991 for a description of these 
simulations). The typical amount of vorticity at the 10 to 15 /i~^Mpc scale for which the rms of the density is below 
0.5 was found to be about 0.2 (in units of Hq). This is well above the values we have estimated in this paper. Though 
it is quite possible that these numerical measurements are spoiled by noise, we do not expect that it could account 
for all the discrepancy between the measured and the predicted vorticities (as suggested by the relative the scatter 
between the two methods suggested in Bernardeau & Van de Weygaert, 1996). 

There are various possible explanations for such discrepancies. It could of course arise from the fact that the 
vorticity at large-scales does not spring from the rare and large caustics but from small scale multi-steaming events 
that cascade towards the larger scales. Such a scenario cannot be excluded but is difficult to investigate by means of 
analytic calculations. It is also possible that the A-body simulations do not address properly the physics of the large 
scales multi-streaming. In particular the two-body interactions should in principle be negligible, a property which 
seems to be hardly satisfied in current A^-body simulations. This shortcoming has been raised by Suisalu & Saar 
(1995), Steinmetz & White (1997) and more specifically by Splinter et al. (1998), where they examine the outcome of 
the planar singularity in phase space. They have found in particular that in classical algorithms the particle's velocity 
dispersions are incorrectly large in all directions. These could turn out to be a major unphysical source of vorticity 
(since the Lagrangian time derivative of the vorticity scales like the curl of the divergence of the velocity anisotropics). 
Specific numerical experiments, that follow for instance the initial density profiles given in this paper, should be carried 
to address this problem more carefully. 

5. Discussion and Conclusions 

We have estimated, within the framework of the gravitational instability scenario, the amount of vorticity generated 
after the first shell crossings in large-scales caustics. The calculations relied on the Zel'dovich approximation which 
yields estimates of the characteristics of the largest caustics and allows explicit calculation of their vorticity content. 
This analysis corresponds to one of the first attempts to investigate the properties of cosmological density perturbations 
beyond first shell-crossing. The previous investigations (Fillmore & Goldreich 1984, Bertschinger 1985) were carried 
out for spherically symmetric systems only, and obviously do not address the physics of vorticity generation. The only 
other means of investigation for this regime is numerical A^-body simulations. 

We found that large scales caustics can provide only an extremely low contribution to the vorticity at scales of 10 
to 15/i^^Mpc. This contribution could be significant only at relatively small scales, when the variance reaches values of 
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Fig. 14. Pr,.{> u>s) in three dimensions for scales char- 
acterized by a a{Rs) of 0.5 (thick lines) and 1 (thin lines) 
and for a n = —2 (solid line), n = —1.5 (long dash line), 
and n = —1 (short dashed line) power spectrum. The 
Fig. 13. Pr^{> los) in two dimensions for scales charac- filled and open circles correspond respectively to the mea- 
terized by a a{Rs) of 0.5 (thick lines) and 1 (thin lines) sured integrated PDF in a CDM simulation at 15/i^^Mpc 
and for a n = —1.5 (solid line), n — ~1 (long dash line), scale with the "Delaunay" or "Voronoi" methods (see 
and n = —0.5 (short dashed line) power spectrum. Bernardeau & Van de Weygaert 1996). 



a few units. This effect is even more important in three dimensions, the difference arising mainly from the coefficient in 
the exponential cut-off. It is therefore unlikely that these caustics can have produced a significant effect on the velocity 
at large scales. In view of these results, it is amply justified to assume that the velocity remains potential down to 
very small scales, i.e. typically the cluster scale at which it is then more natural to expect the multi-streaming regime 
(not only three-flow regime) to play an important role. 

This result provides a complementary view to the picture developed by Doroshkevich (1970) describing the emer- 
gence of galaxy angular momentum from small-scale torque interactions between protogalaxies (a prediction subse- 
quently checked by White (1984), and examined in more detail by Catelan & Theuns, (1996 and 1997)). We rather 
explore here the large scale coherence of the vorticity field that may emerge in a hierarchical scenario from scale much 
larger than the galactic size. The effects we are exploring here does not originate from the two-body interaction of 
haloes as in the picture developed by Doroshkevich, but from the possible existence of large scale coherent vorticity 
field. The conclusion of our work is however that the efficiency with which the large-scale structure caustics generate 
vorticity is rather low. Therefore these results do not really challenge the fact that the small scales interactions should 
indeed be the dominant contribution to the actual galactic angular momenta. 

As a consequence, we do not expect either a significant correlation of the angular momenta at large scale. In 
particular the vorticity field generated in caustics does not seem to be able to induce a significant large scale correlation 
of the galactic shapes which would have been desastruous for weak lensing measurements^. 

Let us reframe this calculation in the context of perturbation theory which has triggered some interest in the last 
few years as a tool to investigate the quasi-linear growth of structures. One key assumption in these techniques is 
that the velocity field is assumed to form a single potential flow. The detailed description of the properties of the first 
singularities is by essence not accessible to this theory: such singularities cannot be "seen" through Taylor expansions 
of the initial fields. In this context it was unclear whether the back reaction of the small scales multi-streaming regime 



^ In these measurements background galaxy shapes are assumed to be totally uncorrelated in the source plane, the observed 
correlation being interpreted as entirely due to gravitational lens effects. 
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on the larger scales (which were thought to be adequately described by perturbation theory) could affect the results on 
those scales. Such effects are partially explored here where we find that the impact of the first multistreaming regions is 
rather low on larger scales. Our results therefore support the idea that the large scales velocity field can be accurately 
described by potential flows and support our views on the validity domain of perturbation theory calculations. 

In the course of this derivation we have made various assumptions. We followed in essence the approach pioneered 
by Press & Schechter (1974) for the mass distribution of virialized objects by trying to identify in the initial density 
field the density fluctuations that contribute mostly to the large-scales vorticity. The calculations have been designed 
to be as accurate as possible in the rare event limit, an approximation which turned out to be crucial at various stages 
of the argument. 

— The above estimation relies heavily on the assumption that the caustics only contribute to large-scale vorticity 
independently of each other. In other words it is assumed that the caustics do not overlap. Moreover the dynamical 
evolution of one caustic is taken to be well-described by the evolution of the caustic having the mean profile. 
This can be approximately true only in the rare event limit since otherwise it is likely that the substructures and 
its environment will change the dynamical evolution of the caustics. Although it is clear that, in the regime we 
investigated, the caustics are rare enough not to overlap, the effects of substructure are more difficult to investigate. 
In particular we have outlined some local features (3.3.1) of the vorticity maps that we think are unlikely to survive 
the existence of substructures. 

— The typical caustics are characterized in this rare event limit. For instance the values of a and b were found to be 
all the more peaked to given values as the corresponding events are rare. We have then estimated the vorticity such 
caustics generate while assuming that slightly different geometries are unlikely to produce very different results. 
This assumption is somewhat suspicious, since it might turn out that slightly different geometries could produce 
more vorticities, and thus change the exact position of the cut-off. We do not expect however that the conclusions 
we have reached could be changed drastically in this manner. 

— The contributions of each caustics to Pj^^ (> a;,,) have also been calculated in the rare event limit. This is in practice 
a very useful approximation on large scales since it is then natural to expect the entire distribution to be dominated 
by a unique value of Amax- 

— We have finally deliberately simplified the spatial distribution of the vorticity within the caustics. Since in the rare 
event limit it is natural to expect that the Lagrangian scales of the caustics are much smaller than the smoothing 
scale this detailed arrangement should be of little relevance. It certainly should not affect the scaling laws as only 
the value of the overall factor Vq will change, and this has little bearing on our conclusions. 

On top of the rare event limit approximation, we have also made a dramatic simplification by using the Zel'dovich 
approximation throughout. This is certainly a secure assumption before the first shell-crossing since the geometries 
that we have investigated were rather sheet-like structures (and the Zel'dovich approximation is exact in ID dynamics). 
After the first shell-crossing however, the back reaction of the large over-densities that are created could possibly affect 
the velocity field. However we do not expect that this effect should be very large so long as Amax is moderately small 
(up to about 1.5), since before then the initial inertial movement should dominate. Later on, matter is expected to 
bounce back to the center of the caustics. Whether the vorticity content is then amplified or reduced remains an open 
question. 



CP wishes to thank J.F. Sygnet, D. Pogosyan, S. Colombi and J.R. Bond for useful conversations. Funding from 
the Swiss NF is gratefully acknowledged. 
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A. Average profile of an a-spherical constrained random field 

A.l. General Formula 

Let us evaluate here the average profile of an a-sphcrical constrained random field in both 2 and 3D. Similar calculations 
as those presented in this Appendix have been investigated by Bardeen et al. (1986) for the 3D field and by Bond 
& Efstathiou (1987) for 2D fields. But, here, instead of the second order derivative of the density field, we consider 
instead the deformation tensor corresponding to second order derivatives of the potential. We also investigate the 
global properties that such constraints induce on the density field. 

Consider a random density field, in either 2D or 3D, having fiuctuations following a Gaussian statistics. It is then 
entirely determined, in a statistical sense, by the shape of its power spectrum, P{k). Recall that P(k) is defined from 
the Fourier transform of the density field, 

(5(k) = /"d3xexp(ik.x) (5(x), with ((5(k)(5(k')) =(5Dirac(k + k')P(fc), (Al) 

where the brackets L') stands for the ensemble average of the random variables. Let us calculate the expectation value 
of (5(k) when a local constraint has been set in order to create an a-spherical perturbation. To set such a constraint, we 
have chosen to consider the deformation tensor of the density field smoothed at a given scale Rl- This tensor reads, 

0,,, = y d'k <5(k) WoikRL) ^. (A2) 

Note that the local smoothed density is given by the trace of this tensor. The chosen window function Wd in Fourier 
space corresponds to a top-hat filter in real space and it reads, 

^2(fc)-2:|^ in 2D, W^{k)^i./^^-^0^ in iD, (A3) 

where J^ are the Bessel functions of index v. The matrix 0i.j is now set to be equal to a given constraint. It is obviously 
possible to choose the axis so that this constraint is a diagonal matrix with eigenvalues {\i),i — 1,D. The elements of 
the matrix 0i.j and (5(k) form a Gaussian random vector, 

Vc = {S{k),(j)i^l,. . . ,(j)D,D,4>l,2,- ■ ■ ,4>1,D, 02,2, • • • ,0_D,_D-l) , (A4) 

and the desired expectation value of S(k) is directly related to the cross-correlation matrix of the components of this 
vector. Consider the matrix Ca_b with a = 0, • • • D{D + l)/2 and 6 = 0, • • • D{D + l)/2, so that 



Co,o = (^(k) <5(k)) = P(fc) , (A5) 

i 3 
"fc2" 



Ca,o = (Sik) 4,,) = Pik) WoikEL) ^ , (A6) 



Ca,b = {cb.^,A'.r)=J^'^Pik)WUkRL) ^'^f:'^'\ (A7) 
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where the indices i,j (respectively i',j') for the matrix elements (j>ij corresponds to the (a+1)"^ (respectively (6+1)*'^) 
component of Vc- For a given spectrum these quantities are easily calculated and are given in the following subsections 
for power law spectrum in resp. 2 and 3 dimensions. The distribution function of the components of the vector Vc then 
reads in terms of Eq. (A7), 



p{Vc) dVc = exp 



lJ2ic'')a,,y^^y'^b 



dVr 



[2^Det(C)]i/2+i3(i3+i)/4- 



(A8) 



The expectation value of 5{h) is given by the ratio 

.expec.., ^ ^ J dS jk) 6 jk) p{V c) 

^ ' Jd5{k)p{Vc) ' 

A straightforward calculation shows that this quantity is given by 
D (C^^ 
,^1 \^ Jo,o 



(A9) 



(AlO) 



Note that the further constraint that the first derivative of the density field should be zero (so that the point xq is 
actually located on a maximum of the density field) would not change the resulting expression of (5°''P°^(k) since the 
cross correlation of the first order derivatives with any other involved quantities identically vanish. 

A. 2. The 2D profile 

In 2 dimensions we have 



(Cofl 


Co.l 


Co, 2 


Co.3\ 


Co.i 


3a^/8 


a^/8 





Co, 2 


ai/8 


3aoV8 





VCo,2 








'^oVs; 



Ca,i 



with the variance of the smoothed density field, ctq, given by 

al^ fd'^kPik) WlikRL). 

The required elements of the inverse of this matrix are given by 
1 



(All) 



(C-) 



0,0 



64 



(^-^)o,i = - 



</Det(C) , 

Co,l Co. 2 Co. 3 



ial, 











afu 



1 



(Co,2 ^ 3Co,i) ctq 



(C-) 



0,2 



Co,l Co^2 Co, 3 

crll8 all 8 



64 Det(C) 64 Det(C) 

1 _ (Co,i-3Co,2)c7^ 



all^ 64Det(C) 
As a result, Eq. (AlO) becomes here 
P{k)W2{kRL) 



64 Det(C) 



(A12) 

(A13) 
(A14) 

(A15) 



^oxpoc.^l^^ 



(Ai + A2 + 2cos(20)[Ai-A2]) , 



(A16) 



where the angle 6 were chosen so that 

ki/k — cos(6') , k2/k = sin(9). 

9 the angle between a given vector and the eigenvector associated to the first eigenvalue. 
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A. 3. The 3D profile 

In 3 dimensions the matrix C reads, 

/3 1 1 \ 

13 1 
1 1 3 

1 
10 
V 00 1/ 

From this expression of the matrix of the eross correlations it is quite straightforward to re-express Eq. (AlO) as 

W{k)W3ikRL) 



C = 



Co,o . . . Cq,q ^ 
: D 


, with 


D = 


1^ 


Co, 6 / 









(A17) 



(5°'=P°=-(k) 



(Ai[k2 + k2 - 4k2]+ A2[k2 + k2 - 4k2] + \^\V\ + k^ - 4k2]) 



When the coordinates of the wave vector are expressed in terms of the angles Bk and (\)u , defined by 
ki — k sm(6k) cos{(f>k) k2 ~ k sni{6k) sai{(l)k) and k^ — k cos{6k) ■ 
Eq. (A18) becomes 

3P{k)W3ikRL) 



^cxpcc.(-j^^ 



(Al + A2 + 6A3) (1 + a cos(26lfe) + b cos(2(/)fe)[l + cos(26lfe)]) 



(A18) 



(A19) 



where a and b are specific combinations of the eigenvalues, 
2A3 — Al — A2 , , ^ Al — A2 



Al + A2 + 6A3 



and b 



Al + A2 + 6A3 



(A20) 



B. The DF of the eigenvalues of the local deformation tensor 

The derivation of the distribution function of the eigenvalues of the local deformation tensor was carried in 3D by 
Doroshkevich (1970). We extend here the calculation to the 2D case (for which the calculations are straightforward). 
Starting with equation (All) - the cross-correlations between the elements of the deformation tensors, one can easily 
get the expression of the joint distribution function of the deformation tensor elements. 



p(0i,i' '^1,2, 02,2) d01,l d01,2 

The change of variables, 
A, 



8 d^i^i d4>i^2 d02,2 



(27r)3/2 



exp 



- 2 (^'^M + ^'^1^2 + 302,2 ~ 201,202,2) 



91.1 



P2,2 



A 



2 2 ^ ' 



"2,2 



r+' 



*?.2, 



2 2 

allows us to introduce the eigenvalues of the matrix. The Jacobian J of this transformation is given by 

, 4>1,1~<P2,2 



J-' = 



2VA 



01,1 — <P2, 2 

2v^ 



1 

As a result we have 

p(A+,A_,0i_2) dA+ dA_ d0i,2 = 



pi, 1—92,2 

2VA 

i>l,l + 4>2,2 

2v^ 





v/l-4 0i,2/A 



BdA , dA_ 



"1.2 



(2^)3/2^3 ^l_4 0i,2/A 



exp 



-4(^JN4J2) 
"0 ^ 



with 

Ji = A+ + A_ , , and J2 = A+ A_ 

The integration over 0i^2 yields 

_ pi dA+ dA_ 



p(A+,A_)dA+dA_ = ./^ 



A, -A- 



exp 



uin-^j. 



Note that if A+ is a priori assumed to be greater than A_ the distribution should be multiplied by 2. 



(Bl) 



(B2) 



(B3) 



(B4) 

(B5) 
(B6) 
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C. Estimation of Pr^{> u)s) 

In this Appendix we estimate the probabihty Pr^{> lUs) that a sphere of radius Rs contains an integrated vorticity 
larger than lUs- In order to account for caustics of all sizes we argued in the main text that Pr^{> uJs) was well 
approximated by 



PRsi>^s 



max 

Rl 



d \ nji^{\) Vca,us.{RL,Rs,{^i}T^s) 



(CI) 



We will now show that the maximum is indeed given by caustics of size of the order of Rs and approximate this 
integral in 2 and 3D. To simplify further Eq. (CI), note first that the distribution function of the eigenvalues is peaked 
in a given geometry (i.e. a ~ 1, and 6 ~ in 3D) for rare caustics (large values of Amax)- Therefore the integral in 
Eq. (CI) will be dominated by caustics of this geometry and the factor Fcaus. can be taken at this point while carrying 
the integration over the other two eigenvalues. As a result we have 



Pr^(> LUs) — max 

Rl 



dA. 



T^maxl^'^D 



/\max 


\ '^o(Amax) K:aus.(^L,-Rs, Amax, Ws) 


^(Rl), 


' RE \ 



(C2) 



This integral runs from 1 to infinity since the caustics exist only when Amax is greater than 1. The evaluation of 
Eq. (C2) requires insights into the function Vcaus.- Although there are no real qualitative changes between the the 2D 
and 3D cases, we now proceed with the computation of Eq. (C2) by distinguishing the two geometries for the sake of 
clarity. 



C.l. The 2D statistics 

Recall that the integral Eq. (CI) will be dominated by the rare even tail, and thus by the lowest value of Amax that 
contributes to the integral. In other words, when considering a given caustic characterized by its Lagrangian scale Rl, 
one should wait long enough so that it has grown sufhciently in order to contribute after sampling a vorticity larger 

is non zero: 



PrA>^s) 



max 

Rl 



Cl^niax Pniaxl^'^max j 



(0) 

a> 

An 



than ujs- For each Rl therefore corresponds Ainax(fiL), the lowest value of Amax for which Vc 

2 



fT(i?L) 



'^o(Amax) Venus. \Rl, Rs, Aniaxj ^s) 



Rl 



(C3) 



The lower bound Amax(^L) is reached as soon as Wquad. is larger than tt i?^ lJs' the largest possible value of the 
integrated vorticity in a cell of a given radius. It is therefore implicitly defined by 



^quad. 
TT Rl 



■ ^M 



im 



RL 

ttRI 



LOo{\^X{^s,RL)-ir. 



(C4) 



Assuming that Vcaus. does not contain any exponential cutoff, and assuming that Amax is in the rare event tail, Eq. (C2) 
can be approximatively re-expressed as 



Pr, (> LjJs 



max 

Rl 



0.56 



A 



(0) 



(y{RL) 



exp 




no Vca.us.{Rs,RL,^ 



(0) 



,^s 



Rl 



(C5) 



when using Eq. (11) for the distribution function of Amax, integrating by part and dropping the residual integral for 
large enough Amax/c(-RL) (see Appendix E for details). This maximum with respect to Rl is then approximated by the 
minimum of the argument of the exponential, Amax(^L)/o'(i?L)7 where the minimum in the facto taken with respect 
to Amax since cr(i?L) can be thought of a function of Amax via Eqs. (22) and (C4). This minimum can de facto be 
expressed independently of Rs . It reads 



A(") 



4-a(n + 2) 



(C6) 



Once Amax is fixed the geometry of the caustic which will contribute most to Pr^{> uJs) is entirely specified. The 
condition for the existence of a minimum defining Amax is that a{n + 2) < 4, and it is satisfied for all considered 
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cases (see Table (1)). This implies that we are investigating a regime where the integral Eq. (C2) is not dominated 
by arbitrarily rare caustics - which would have been catastrophic given the assumptions (note that when n is too 
large Amax tend to be quite large thus challenging the validity of quantitative results based upon the Zcl'dovich 
approximation). The resulting value of Rl is 

«, = «, ,/^: fi^i^i^V" . /, «, (^\"\ ,C7) 



The scale factor /„ is given in Table (3) for an Einstein-de Sitter universe (/(fi) ~ 1) and different values of n. 
Completing the calculation of Pr^{> (jJs) involves relating the shape and size of the caustic for the adopted value of 
Amax- These values are derived from the fits (Eq. (47)) and are given in Table (3). Fig. (12) gives Kaus., in units of 
the square of Rl, as a function of the smoothing radius Rg. From Fig. (12) it is easy to see that 

Kaus. ^ Vo Rs Rl , (C8) 

for any values and n; the corresponding values of Vq are given in Table 3. Putting Eq. (C8) into Eq. (C3), using 
Eqs. (C6), (C7) yields for the sought distribution 




Pr^ (> a;,) ^ 0.56 no Vo I ^^^ I f:+' uj^-^^^^'^ exp 




«+2^(n+2)/2 



(C9) 



Note that the power of ujg in the exponential is rather weak. The cut-off is nonetheless strong in the regime of interest 
because of the leading coefficient. Equation (C9) is illustrated on Fig. (13) and discussed in the main text. 



Table 3. Parameters of interest for the 2D caustics: the power index, n, the critical time Amix, the radial extension e '% depth 
d}^' in units of Rl , scale factor fs as well as the values of no and Vo for the critical caustics. 



n 


),(0) 


d(") 


e(0) 


/. 


no 


Vo 


-1.5 


1.31 


0.17 


1.34 


0.30 


0.018 


0.9 


-1 


1.67 


0.40 


1.33 


0.95 


0.023 


1.8 


-0.5 


2.15 


0.90 


1.36 


1.25 


0.009 


3.4 



Table 4. Parameters of interest for the 3D caustics; the power index, n, the critical times A^ax, the scale factor /^ in the two 
regimes {Rs < e/3 in parentheses) with radial extension e'"', depth d'"' in units of Rl as well as the values of no and Vo that 
enter the final expressions. 



n 


\ + 

^max 


(Amax) 


/.+ 


Us) 


d(») 


e(") 


no 


Vb 


-2. 


1.41 


(1.47) 


2.46 


(2.09) 


0.18 


1.04 


0.18 


0.96 


-1.5 


1.63 


(1.79) 


2.10 


(1.58) 


0.28 


1.05 


0.14 


1.84 


-1. 


1.84 


(2.15) 


1.78 


(1.17) 


0.42 


1.07 


0.064 


3.18 



C.2. the 3D statistics 

The threshold on Amax, from which the caustics start to contribute at a given scale Rs depends on the adopted 
description for the local vorticity. We assume here as mentioned in Sect. 4.3.2 that the total vorticity is localized on 
two rings of radius e/3 each, distant of 2 d/i of each other. They are assumed to bear opposite lineic (and uniformly 
distributed) vorticitics; in order to get a consistent answer for the integrated vorticity in a quadrant, we should have 

o.,„. = ^^:^3ii^. (CIO) 
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^J^L Fig. 16. The loci of the centres of spheres contributing lOs 

in the range [0)^(1 — 6^/2), o)^[. The dashed arrow points 
to a centre of such a sphere, and defines the running angle. 
Fig. 15. The function Kaus. , in units of the square of Rl , 0, mentioned in Eq. (D2). The two (cylindrically symmet- 
as a function of the smoothing radius in 2D. The solid ric) shaded regions correspond to the loci of the centre of 
line corresponds to the case n — —1.5, the dashed line to spheres capturing almost half a ring and all or none of 
n — —1 and the long dashed line to n = —0.5. In all cases the other. Two examples of such spheres are displayed for 
the geometry of the caustic is fixed by Amax ~ Amax. either case. 

The maximum vorticity that can be encompassed in a sphere then depends on its radius Rs- If Rs is larger than the 
radius of the rings e/3, it is possible to have half of a ring in a sphere (while the other ring does not intersect it at 
all), so that the values of A^ax (for which the maximum vorticity is sampled) is given by 



^^ = ^^4^^-M = |^-o(A„.a.-ir/(f^)^, if Rs>e/3. 



(Cll) 



If on the other hand R^ is smaller than e/3 then only a fraction of the half ring can be put in the sphere and we have 
instead 



c^, = 2i?,o;H„. 4^^ ^M = |^^ (Amax ~ir-"V(r!)^, if i?. < e/3. 



(C12) 



Now the local behaviour of T^aus. near its takeoff value is well represented (as argued below and demonstrated in 
Appendix D for large enough Rs) as a function of A^ax by 



Vc,us.iRL,Rs, A„,ax, i^s)^ fo K (C, Rl,Rs, Amax) - t^s] d^C ~ i^L^^'K) (Amax - ^^LV : 



(C13) 



Using Eq. (17) and (52) for the distribution function Pmax(Amax)7 changing integration variable from u = Amax/c to 
Amax + u/Amax and dropping the residual integral for large enough Amax/o'(i?L) (see Appendix E for details) yields for 
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Eq. (C3) : 
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From Eq. (22) and (Cll), (C12), the minimum of the argument of the exponential corresponds to: 
6 .„„,..._ 4 



A 



6 — a{n + 3) 



if Rs > e/3, and Aj^ 



4 — (a — ae){n + 3) 



if Rs < e/3, 



(C14) 



(C15) 



which assumes that a(n + 3) < 6 (resp. {a — ae)(n + 3) < 4 ), both conditions being satisfied for all values of n 
considered. The corresponding scaling relations between Rj^ and Rs are given by 



Rh — ,fs Rs 



\fm) 



1/3 



/ \ 1/2 

if Rs > e/3 , or i?L - ./7 Rs ( -^ j H Rs < e/3 . 



(C16) 



The scale factors ff - derived from the fits (Eq. (47)) - are given in Table (4) for an Einstein-de Sitter universe 
(/(ri) = 1) and different values of n. Interestingly, as long as uJs is not too large the condition Rl > e/3 is always 
satisfied. In practice at scales of about 10 to 15/i~^Mpc the measured vorticity tUs is expected to be indeed at most 
of a few tenth (Bernardeau & van de Weygaert, 1995). It is therefore always fair to assume that we are in the regime 
where Rs > e/3 which is the regime investigated hereafter. 

Completing the calculation of Pr^{> uJs) requires evaluating the corresponding no, 7 and Vq. The value of rig is 
entirely determined by the geometry of the caustics and is given in the Tables 1 and 2. The behaviour of V^aus. a-s it 
departs from zero as a function of lOs for the critical ratios oi Rs, e and d is locally well fitted as a function of ujs by a 
power law of the form 

(C17) 



K;aus.(Amax,^s) — Uq Rl Rs { ^ ^ 



t^M(-^max) 



where w^f is the threshold value of ojs (Eq. (Cll)). This expression is vahd when iVs is close to its threshold value. On 

it is possible to relate the variation of A^ax to the variations of w^. We can then rewrite 



^M' 



the critical line, 

Eq. (C17) as a function of the difference between Amax and the critical value Amax, assuming this departure is small, 



Kaus. (Amax, Ws) ~ i?L -RgVo (Amax - A^],^)'^ i with Vq 



(A 



(0) 



1)7 



(C18) 



Since Rl/Rs is only a function of n and ujs, so are Vq and 7. In practice we take the asymptotic values of Vq and 
7 given in Appendix D and corresponding to the limit Rg <C Rl- Putting Eq. (C18) into Eq. (C14), using Eq. (C15) 
- (C17) and (D4) yields for the vorticity distribution 



PB,A>^s)^OA8noVo 



«(0) \V2 

-^max \ 



CT(i?. 



(13 + 7n) (13 + 7n 

4 , , 12 



exp 



^max 

(t(Rs) 



f"+3^(«+3)/3 



(C19) 



Equation (C19) is illustrated on Fig. (14) and discussed in the main text. 



D. Asymptotic behaviour of Vcaust. in 3D 

For large enough Rg we derive here an asymptotic analytic expression for Vcaust.- Let us first estimate geometrically 
the volume in space contributing almost wj^ to Vcaust.- The corresponding contribution is the sum of two volumes 
given by the shaded area in Fig. (16), corresponding to the loci of the centers of spheres which capture almost half a 
ring and not the other, or which capture completely one ring and almost half of the other. In the asymptotic limit, as 
e/ Rs -^ 0, the element of volume is an infinitely thin strip and both contributions become equal since 6 -^ —6'. The 
area corresponding to these loci can be evaluated algebraically as follow: let us call e the projected ring segment by 
which a sampling sphere of radius Rg fails to encompass a ring diameter 2e/3; it follows that the ratio of oJs to w^, 
is given by 






(1 



(Dl) 



M 
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On the other hand, for a given direction for the sphere centre given by cos{9) — ^, within the solid angle 27rd/i, the 
volume clement (encompassed by the two shifted spheres capturing ojg in the range [^'-'^^(l — £^/2), ^m\) is given by 

2e 



47r Rlesm^ OAO = Sn-R^.e^/l - m^ d^ . 
o / o 



Summing over all possible directions {i.e. before intersecting the second ring) yields 



(D2) 



g r g(0) 

87r-i?2e / y/i - ^2 d^ = Sn—-RLR'^^eJ, where 



fJ-a = 



Mo 






21 -1/2 



(D3) 



Accounting for the summation over the two configurations (half a ring captured or a full + one half ring captured) , 
using Eq. (Dl) to eliminate e, we finally get for large enough Rs 



^Caust. = 16V27ri?i Rl^ ( 1 - ^ 

■i \ Ujjj 



1/2 



1 r- e'^^' 

therefore 700 = :r and U^ = 16v27r— - J . 



(D4) 



E. Rare event approximation 

Consider an integral of the form 



x^ (x — a)'^ exp(— 6x ) dx . 



(El) 



Changing variable to x — a + u/{2ab) Eq. (El) reads 

00 

U \/3 



i.5L,.„_,„^,/(^)\. 



1 



2a% 



exp 



2ba' 



exp(— u) du . 



(E2) 



For large enough a the square brace in Eq. (E2) is well approximated by 1 yielding for Eq. (E2) 
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j-r(7 + l) exp(-fea" 
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Eq. (C5) is a special case of Eq. (El) with x = Amax/cr, a = Ajnax, 7 = 0, /3 = 3 and b = 4/3, while Eq. (C14) 
corresponds to /? = 5, and b = 5/2. Note that the 7 = approximant can be deduced directly by integration by parts. 
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